Investigating Second Eigenvalue

  if (require("PageRank")) {
      library(PageRank)
    }else{
      devtools::install_github("ryangreenup/PageRank")
      library(PageRank)
    }
Loading required package: PageRank
  library(pacman)
  pacman::p_load(PageRank, devtools, Matrix, igraph, mise, tidyverse, rgl, latex2exp)
#  mise()

Looking at Density

Constants

Define some constants

p    <- seq(from = 0.01, to = 0.99, length.out = 10)
beta <- seq(from = 1, to = 20, length.out = 40)
sz <- seq(from = 100, to = 10, length.out = 5)
input_var <- expand.grid("p" = p, "beta" = beta, "size" = sz)
input_var

Function to Build Graph

random_graph <- function(p, beta, size) {
      g1 <- igraph::erdos.renyi.game(n = size, p)
      A <- igraph::get.adjacency(g1) # Row to column
      A <- Matrix::t(A)

      A_dens <- mean(A)
      T      <- PageRank::power_walk_prob_trans(A, beta = beta)
      tr     <- sum(diag(T))
      e2     <- eigen(T, only.values = TRUE)$values[2] # R orders by descending magnitude
      return(c(abs(e2), mean(A), tr))
}

Return results

Map the function

nc <- length(random_graph(1, 1, 1))
Y <- matrix(ncol = nc, nrow = nrow(input_var))
for (i in 1:nrow(input_var)) {
  X <- as.vector(input_var[i,])
  Y[i,] <-  random_graph(X$p, X$beta, X$size)
  print(i/nrow(input_var))
}
[1] 5e-04
[1] 0.001
[1] 0.0015
[1] 0.002
[1] 0.0025
[1] 0.003
[1] 0.0035
[1] 0.004
[1] 0.0045

Plot Results

pairs(data)

cor(data)
                        p          beta          size eigenvalue2
p            1.000000e+00 -9.147063e-21  0.000000e+00  -0.5312247
beta        -9.147063e-21  1.000000e+00  8.159342e-21   0.4406453
size         0.000000e+00  8.159342e-21  1.000000e+00  -0.4308528
eigenvalue2 -5.312247e-01  4.406453e-01 -4.308528e-01   1.0000000
A_dens       9.999629e-01 -1.040997e-04  2.837892e-03  -0.5324257
trace       -5.318629e-01 -5.943376e-01 -2.002601e-03  -0.1202697
                   A_dens        trace
p            0.9999628681 -0.531862883
beta        -0.0001040997 -0.594337559
size         0.0028378923 -0.002002601
eigenvalue2 -0.5324256616 -0.120269676
A_dens       1.0000000000 -0.531765466
trace       -0.5317654658  1.000000000
library(corrplot)
corrplot 0.84 loaded
cormat = cor(data, method = 'spearman')
corrplot(cormat, method = "ellipse", type = "lower")

names(data)
[1] "p"           "beta"        "size"        "eigenvalue2"
[5] "A_dens"      "trace"      

Let’s look at a 3d Output

# plot3d(data$beta, data$A_dens, data$eigenvalue2, size = 2)

names(data)
[1] "p"           "beta"        "size"        "eigenvalue2"
[5] "A_dens"      "trace"      
library(plotly)

d <- data[sample(1:nrow(data), 1000),]
d <- data
# d$beta <- log(d$beta)

fig <- plot_ly(d, x = ~A_dens, y = ~beta, z = ~eigenvalue2)
fig <- fig %>% add_markers(size = 1)
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Density'),
                     yaxis = list(title = 'Beta'),
                     zaxis = list(title = 'E2')))

fig

NA

names(data)
[1] "p"           "beta"        "size"        "eigenvalue2"
[5] "A_dens"      "trace"      
library(plotly)

d <- data[sample(1:nrow(data), 1000),]

fig <- plot_ly(d, x = ~A_dens, y = ~trace, z = ~eigenvalue2)
fig <- fig %>% add_markers(size = 1)
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Density'),
                     yaxis = list(title = 'trace'),
                     zaxis = list(title = 'E2')))

fig

NA

Clearly I should be able to model this.

Let’s look at density, it’s continuous so I should be able to take splices of A_dens.

A_dens is dependent on p, so I’ll just use different ps

p    <- seq(from = 0.1, to = 0.95, length.out = 5)
beta <- seq(from = 1, to = 10, length.out = 1000)
size = 100
input_var <- expand.grid("p" = p, "beta" = beta, "size" = size)
input_var

nc <- length(random_graph(1, 1, 1))
Y <- matrix(ncol = nc, nrow = nrow(input_var))
for (i in 1:nrow(input_var)) {
  X <- as.vector(input_var[i,])
  Y[i,] <-  random_graph(X$p, X$beta, X$size)
  print(i/nrow(input_var))
}
if (sum(abs(Y) != abs(Re(Y))) == 0) {
  Y <- Re(Y)
}
nrow(input_var)
nrow(Y)
Y <- as.data.frame(Y); colnames(Y) <- c("eigenvalue2", "trace")
(data2 <- cbind(input_var, Y)) %>% head()
data2$p <- factor(data2$p)
ggplot(data2[30:nrow(data2),], mapping = aes(col = p, x = beta, y = eigenvalue2)) +
  geom_point(size = 0.5) +
  stat_smooth() +
  scale_size_continuous(range = c(0.1,1)) +
  labs(x = "Beta", y = TeX("Second Eigenvalue"), title = TeX("Second Eigenvalue given Matrix Density") ) +
  guides(col = guide_legend("Link Density"))  +
  theme_bw()

These look logarithmic, let’s do a log transform:


ggplot(data2[30:nrow(data2),], mapping = aes(col = p, x = log(beta), y = eigenvalue2)) +
  geom_point(size = 0.5) +
  stat_smooth() +
  scale_size_continuous(range = c(0.1,1)) +
  labs(x = "Beta", y = TeX("Second Eigenvalue"), title = TeX("Second Eigenvalue given Matrix Density") ) +
  guides(col = guide_legend("Link Density"))  +
  theme_bw()

Link density is a function of size so E2ln(Beta) is still pretty useful.

it seems that this is very nearly linear for low densities, the web would have a low density, for example ba graphs appear bounded above by 1/2:

n_vec <- 0:20
m <- 3
y <- c()
for (n in n_vec) {
    x <- sample_pa(n = n, power = 3, m = m) %>% 
      get.adjacency() %>% 
      mean()
    y <- c(y, x)
}
plot(n_vec, y)
length(n_vec)
length(y)

So considering small Densities is wise, let’s have a look.

So at this stage I would like to fit a model over what we have in the first plotly diagram, but, I’m going to look at ba graphs in case there is something more insightful.


ggplot(data2[30:nrow(data2),], mapping = aes(col = factor(p), x = log(beta*p), y = eigenvalue2)) +
  geom_point(size = 0.5) +
  stat_smooth() +
  scale_size_continuous(range = c(0.1,1)) +
  labs(x = "Beta", y = TeX("Second Eigenvalue"), title = TeX("Second Eigenvalue given Matrix Density") ) +
  
  guides(col = guide_legend("Link Density"))  +
  theme_bw()

It seems for very low densities that it is linear and then logarithmic for moderate densities. BA models have density~=m/n so for the BA model we’ll just use beta.

LS0tCnRpdGxlOiAiSW52ZXN0aWdhdGluZyBTZWNvbmQgRWlnZW52YWx1ZSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyBJbnZlc3RpZ2F0aW5nIFNlY29uZCBFaWdlbnZhbHVlCgpgYGB7cn0KICBpZiAocmVxdWlyZSgiUGFnZVJhbmsiKSkgewogICAgICBsaWJyYXJ5KFBhZ2VSYW5rKQogICAgfWVsc2V7CiAgICAgIGRldnRvb2xzOjppbnN0YWxsX2dpdGh1YigicnlhbmdyZWVudXAvUGFnZVJhbmsiKQogICAgICBsaWJyYXJ5KFBhZ2VSYW5rKQogICAgfQoKICBsaWJyYXJ5KHBhY21hbikKICBwYWNtYW46OnBfbG9hZChQYWdlUmFuaywgZGV2dG9vbHMsIE1hdHJpeCwgaWdyYXBoLCBtaXNlLCB0aWR5dmVyc2UsIHJnbCwgbGF0ZXgyZXhwKQojICBtaXNlKCkKYGBgCgoKCiMjIExvb2tpbmcgYXQgRGVuc2l0eQoKIyMjIENvbnN0YW50cwoKRGVmaW5lIHNvbWUgY29uc3RhbnRzCgpgYGB7cn0KcCAgICA8LSBzZXEoZnJvbSA9IDAuMDEsIHRvID0gMC45OSwgbGVuZ3RoLm91dCA9IDEwKQpiZXRhIDwtIHNlcShmcm9tID0gMSwgdG8gPSAyMCwgbGVuZ3RoLm91dCA9IDQpCnN6IDwtIHNlcShmcm9tID0gMTAwLCB0byA9IDEwLCBsZW5ndGgub3V0ID0gNSkKaW5wdXRfdmFyIDwtIGV4cGFuZC5ncmlkKCJwIiA9IHAsICJiZXRhIiA9IGJldGEsICJzaXplIiA9IHN6KQppbnB1dF92YXIKYGBgCgojIyMgRnVuY3Rpb24gdG8gQnVpbGQgR3JhcGgKCmBgYHtyfQpyYW5kb21fZ3JhcGggPC0gZnVuY3Rpb24ocCwgYmV0YSwgc2l6ZSkgewogICAgICBnMSA8LSBpZ3JhcGg6OmVyZG9zLnJlbnlpLmdhbWUobiA9IHNpemUsIHApCiAgICAgIEEgPC0gaWdyYXBoOjpnZXQuYWRqYWNlbmN5KGcxKSAjIFJvdyB0byBjb2x1bW4KICAgICAgQSA8LSBNYXRyaXg6OnQoQSkKCiAgICAgIEFfZGVucyA8LSBtZWFuKEEpCiAgICAgIFQgICAgICA8LSBQYWdlUmFuazo6cG93ZXJfd2Fsa19wcm9iX3RyYW5zKEEsIGJldGEgPSBiZXRhKQogICAgICB0ciAgICAgPC0gc3VtKGRpYWcoVCkpCiAgICAgIGUyICAgICA8LSBlaWdlbihULCBvbmx5LnZhbHVlcyA9IFRSVUUpJHZhbHVlc1syXSAjIFIgb3JkZXJzIGJ5IGRlc2NlbmRpbmcgbWFnbml0dWRlCiAgICAgIHJldHVybihjKGFicyhlMiksIG1lYW4oQSksIHRyKSkKfQpgYGAKCiMjIyBSZXR1cm4gcmVzdWx0cwoKTWFwIHRoZSBmdW5jdGlvbgoKYGBge3J9Cm5jIDwtIGxlbmd0aChyYW5kb21fZ3JhcGgoMSwgMSwgMSkpClkgPC0gbWF0cml4KG5jb2wgPSBuYywgbnJvdyA9IG5yb3coaW5wdXRfdmFyKSkKZm9yIChpIGluIDE6bnJvdyhpbnB1dF92YXIpKSB7CiAgWCA8LSBhcy52ZWN0b3IoaW5wdXRfdmFyW2ksXSkKICBZW2ksXSA8LSAgcmFuZG9tX2dyYXBoKFgkcCwgWCRiZXRhLCBYJHNpemUpCiAgcHJpbnQoaS9ucm93KGlucHV0X3ZhcikpCn0KaWYgKHN1bShhYnMoWSkgIT0gYWJzKFJlKFkpKSkgPT0gMCkgewogIFkgPC0gUmUoWSkKfQpucm93KGlucHV0X3ZhcikKbnJvdyhZKQpZIDwtIGFzLmRhdGEuZnJhbWUoWSk7IGNvbG5hbWVzKFkpIDwtIGMoImVpZ2VudmFsdWUyIiwgIkFfZGVucyIsICJ0cmFjZSIpCihkYXRhIDwtIGNiaW5kKGlucHV0X3ZhciwgWSkpICU+JSBoZWFkKCkKYGBgCgoKIyMjIFBsb3QgUmVzdWx0cwoKCmBgYHtyfQpwYWlycyhkYXRhKQpjb3IoZGF0YSkKbGlicmFyeShjb3JycGxvdCkKY29ybWF0ID0gY29yKGRhdGEsIG1ldGhvZCA9ICdzcGVhcm1hbicpCmNvcnJwbG90KGNvcm1hdCwgbWV0aG9kID0gImVsbGlwc2UiLCB0eXBlID0gImxvd2VyIikKbmFtZXMoZGF0YSkKYGBgCkxldCdzIGxvb2sgYXQgYSAzZCBPdXRwdXQKCmBgYHtyfQojIHBsb3QzZChkYXRhJGJldGEsIGRhdGEkQV9kZW5zLCBkYXRhJGVpZ2VudmFsdWUyLCBzaXplID0gMikKYGBgCgoKCmBgYHtyfQoKbmFtZXMoZGF0YSkKbGlicmFyeShwbG90bHkpCgpkIDwtIGRhdGFbc2FtcGxlKDE6bnJvdyhkYXRhKSwgMTAwMCksXQpkIDwtIGRhdGEKIyBkJGJldGEgPC0gbG9nKGQkYmV0YSkKCmZpZyA8LSBwbG90X2x5KGQsIHggPSB+QV9kZW5zLCB5ID0gfmJldGEsIHogPSB+ZWlnZW52YWx1ZTIpCmZpZyA8LSBmaWcgJT4lIGFkZF9tYXJrZXJzKHNpemUgPSAxKQpmaWcgPC0gZmlnICU+JSBsYXlvdXQoc2NlbmUgPSBsaXN0KHhheGlzID0gbGlzdCh0aXRsZSA9ICdEZW5zaXR5JyksCiAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICdCZXRhJyksCiAgICAgICAgICAgICAgICAgICAgIHpheGlzID0gbGlzdCh0aXRsZSA9ICdFMicpKSkKCmZpZwoKYGBgCgpgYGB7cn0KCm5hbWVzKGRhdGEpCmxpYnJhcnkocGxvdGx5KQoKZCA8LSBkYXRhW3NhbXBsZSgxOm5yb3coZGF0YSksIDEwMDApLF0KCmZpZyA8LSBwbG90X2x5KGQsIHggPSB+QV9kZW5zLCB5ID0gfnRyYWNlLCB6ID0gfmVpZ2VudmFsdWUyKQpmaWcgPC0gZmlnICU+JSBhZGRfbWFya2VycyhzaXplID0gMSkKZmlnIDwtIGZpZyAlPiUgbGF5b3V0KHNjZW5lID0gbGlzdCh4YXhpcyA9IGxpc3QodGl0bGUgPSAnRGVuc2l0eScpLAogICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAndHJhY2UnKSwKICAgICAgICAgICAgICAgICAgICAgemF4aXMgPSBsaXN0KHRpdGxlID0gJ0UyJykpKQoKZmlnCgpgYGAKCgpDbGVhcmx5IEkgc2hvdWxkIGJlIGFibGUgdG8gbW9kZWwgdGhpcy4KCkxldCdzIGxvb2sgYXQgZGVuc2l0eSwgaXQncyBjb250aW51b3VzIHNvIEkgc2hvdWxkIGJlIGFibGUgdG8gIHRha2Ugc3BsaWNlcyBvZiBBX2RlbnMuCgpBX2RlbnMgaXMgZGVwZW5kZW50IG9uIHAsIHNvIEknbGwganVzdCB1c2UgZGlmZmVyZW50IHBzCgoKYGBge3J9CnAgICAgPC0gc2VxKGZyb20gPSAwLjEsIHRvID0gMC45NSwgbGVuZ3RoLm91dCA9IDUpCmJldGEgPC0gc2VxKGZyb20gPSAxLCB0byA9IDEwLCBsZW5ndGgub3V0ID0gMTAwMCkKc2l6ZSA9IDEwMAppbnB1dF92YXIgPC0gZXhwYW5kLmdyaWQoInAiID0gcCwgImJldGEiID0gYmV0YSwgInNpemUiID0gc2l6ZSkKaW5wdXRfdmFyCgpuYyA8LSBsZW5ndGgocmFuZG9tX2dyYXBoKDEsIDEsIDEpKQpZIDwtIG1hdHJpeChuY29sID0gbmMsIG5yb3cgPSBucm93KGlucHV0X3ZhcikpCmZvciAoaSBpbiAxOm5yb3coaW5wdXRfdmFyKSkgewogIFggPC0gYXMudmVjdG9yKGlucHV0X3ZhcltpLF0pCiAgWVtpLF0gPC0gIHJhbmRvbV9ncmFwaChYJHAsIFgkYmV0YSwgWCRzaXplKQogIHByaW50KGkvbnJvdyhpbnB1dF92YXIpKQp9CmlmIChzdW0oYWJzKFkpICE9IGFicyhSZShZKSkpID09IDApIHsKICBZIDwtIFJlKFkpCn0KbnJvdyhpbnB1dF92YXIpCm5yb3coWSkKWSA8LSBhcy5kYXRhLmZyYW1lKFkpOyBjb2xuYW1lcyhZKSA8LSBjKCJlaWdlbnZhbHVlMiIsICJ0cmFjZSIpCihkYXRhMiA8LSBjYmluZChpbnB1dF92YXIsIFkpKSAlPiUgaGVhZCgpCmRhdGEyJHAgPC0gZmFjdG9yKGRhdGEyJHApCgpgYGAKCgpgYGB7cn0KZ2dwbG90KGRhdGEyWzMwOm5yb3coZGF0YTIpLF0sIG1hcHBpbmcgPSBhZXMoY29sID0gcCwgeCA9IGJldGEsIHkgPSBlaWdlbnZhbHVlMikpICsKICBnZW9tX3BvaW50KHNpemUgPSAwLjUpICsKICBzdGF0X3Ntb290aCgpICsKICBzY2FsZV9zaXplX2NvbnRpbnVvdXMocmFuZ2UgPSBjKDAuMSwxKSkgKwogIGxhYnMoeCA9ICJCZXRhIiwgeSA9IFRlWCgiU2Vjb25kIEVpZ2VudmFsdWUiKSwgdGl0bGUgPSBUZVgoIlNlY29uZCBFaWdlbnZhbHVlIGdpdmVuIE1hdHJpeCBEZW5zaXR5IikgKSArCiAgZ3VpZGVzKGNvbCA9IGd1aWRlX2xlZ2VuZCgiTGluayBEZW5zaXR5IikpICArCiAgdGhlbWVfYncoKQpgYGAKClRoZXNlIGxvb2sgbG9nYXJpdGhtaWMsIGxldCdzIGRvIGEgbG9nIHRyYW5zZm9ybToKCmBgYHtyfQoKZ2dwbG90KGRhdGEyWzMwOm5yb3coZGF0YTIpLF0sIG1hcHBpbmcgPSBhZXMoY29sID0gcCwgeCA9IGxvZyhiZXRhKSwgeSA9IGVpZ2VudmFsdWUyKSkgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDAuNSkgKwogIHN0YXRfc21vb3RoKCkgKwogIHNjYWxlX3NpemVfY29udGludW91cyhyYW5nZSA9IGMoMC4xLDEpKSArCiAgbGFicyh4ID0gIkJldGEiLCB5ID0gVGVYKCJTZWNvbmQgRWlnZW52YWx1ZSIpLCB0aXRsZSA9IFRlWCgiU2Vjb25kIEVpZ2VudmFsdWUgZ2l2ZW4gTWF0cml4IERlbnNpdHkiKSApICsKICBndWlkZXMoY29sID0gZ3VpZGVfbGVnZW5kKCJMaW5rIERlbnNpdHkiKSkgICsKICB0aGVtZV9idygpCmBgYApMaW5rIGRlbnNpdHkgaXMgYSBmdW5jdGlvbiBvZiBzaXplIHNvIEUyXHByb3B0byBsbihCZXRhKSBpcyBzdGlsbCBwcmV0dHkgdXNlZnVsLgoKaXQgc2VlbXMgdGhhdCB0aGlzIGlzIHZlcnkgbmVhcmx5IGxpbmVhciBmb3IgbG93IGRlbnNpdGllcywgdGhlIHdlYiB3b3VsZCBoYXZlIGEgbG93IGRlbnNpdHksIGZvciBleGFtcGxlIGJhIGdyYXBocyBhcHBlYXIgYm91bmRlZCBhYm92ZSBieSAxLzI6CgpgYGB7cn0Kbl92ZWMgPC0gMDoyMAptIDwtIDMKeSA8LSBjKCkKZm9yIChuIGluIG5fdmVjKSB7CiAgICB4IDwtIHNhbXBsZV9wYShuID0gbiwgcG93ZXIgPSAzLCBtID0gbSkgJT4lIAogICAgICBnZXQuYWRqYWNlbmN5KCkgJT4lIAogICAgICBtZWFuKCkKICAgIHkgPC0gYyh5LCB4KQp9CnBsb3Qobl92ZWMsIHkpCmxlbmd0aChuX3ZlYykKbGVuZ3RoKHkpCmBgYAoKU28gY29uc2lkZXJpbmcgc21hbGwgRGVuc2l0aWVzIGlzIHdpc2UsIGxldCdzIGhhdmUgYSBsb29rLgoKU28gYXQgdGhpcyBzdGFnZSBJIHdvdWxkIGxpa2UgdG8gZml0IGEgbW9kZWwgb3ZlciB3aGF0IHdlIGhhdmUgaW4gdGhlIGZpcnN0IHBsb3RseSBkaWFncmFtLCBidXQsIEknbSBnb2luZyB0byBsb29rIGF0IGJhIGdyYXBocyBpbiBjYXNlIHRoZXJlIGlzIHNvbWV0aGluZyBtb3JlIGluc2lnaHRmdWwuCgpgYGB7cn0KCmdncGxvdChkYXRhMlszMDpucm93KGRhdGEyKSxdLCBtYXBwaW5nID0gYWVzKGNvbCA9IGZhY3RvcihwKSwgeCA9IGxvZyhiZXRhKnApLCB5ID0gZWlnZW52YWx1ZTIpKSArCiAgZ2VvbV9wb2ludChzaXplID0gMC41KSArCiAgc3RhdF9zbW9vdGgoKSArCiAgc2NhbGVfc2l6ZV9jb250aW51b3VzKHJhbmdlID0gYygwLjEsMSkpICsKICBsYWJzKHggPSAiQmV0YSIsIHkgPSBUZVgoIlNlY29uZCBFaWdlbnZhbHVlIiksIHRpdGxlID0gVGVYKCJTZWNvbmQgRWlnZW52YWx1ZSBnaXZlbiBNYXRyaXggRGVuc2l0eSIpICkgKwogIAogIGd1aWRlcyhjb2wgPSBndWlkZV9sZWdlbmQoIkxpbmsgRGVuc2l0eSIpKSAgKwogIHRoZW1lX2J3KCkKYGBgCgoKCkl0IHNlZW1zIGZvciB2ZXJ5IGxvdyBkZW5zaXRpZXMgdGhhdCBpdCBpcyBsaW5lYXIgYW5kIHRoZW4gbG9nYXJpdGhtaWMgZm9yCm1vZGVyYXRlIGRlbnNpdGllcy4gQkEgbW9kZWxzIGhhdmUgZGVuc2l0eX49bS9uIHNvIGZvciB0aGUgQkEgbW9kZWwgd2UnbGwKanVzdCB1c2UgYmV0YS4=